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ABSTRACT 


A mixed formulation for calculating static equilibrium and stability 
eigenvalues of nonuniform rotor blades in hover is presented. The static 
equilibrium equations are nonlinear and are solved by an accurate and effi- 
cient collocation method. The linearized perturbation equations are solved 
by a one-step, second-order integration scheme. The numerical results 
correlate very well with published results from a nearly identical stabil- 
ity analysis based on a displacement formulation. Slight differences in 
the results are traced to terms in the equations that relate moments to 
derivatives of rotations. With the present ordering scheme, in which terms 
of the order of squares of rotations are neglected with respect to unity, 
it is not possible to achieve completely equivalent models based on mixed 
and displacement formulations. A study of the one-step methods reveals that 
a second-order Taylor expansion is necessary to achieve good convergence 
for nonuniform rotating blades. Numerical results for a hypothetical non- 
uniform blade. Including the nonlinear static equilibrium solution, were 
obtained with no more effort or computer time than that required for a 
uniform blade with the present analysis. 

1. Introduction 




It has been found that nonlinearities in rotor-blade equations affect 
blade stability [1-4] — especially stability of coupled flap, lead-lag, and 
torsion degrees of freedom [1]. In these references it was found that 
essential nonlinear effects could be retained by perturbing the nonlinear 
equations of motion about the static equilibrium condition and solving for 
the eigenvalues of the linearized perturbation equations. Coefficients of 
the linearized perturbation equations then depend on the solution of the 
nonlinear static equilibrium equations. 

This paper presents a method for stability analysis of nonuniform 
rotating blades with aerodynamic loading that utilizes the solution for the 
nonlinear static equilibrium equations in the linearized eigenvalue problem. 
Methods described in the literature have been limited to solution of 
various restricted versions of this problem [4-9]. These include a modal 
approach [4], an integrating matrix method [5], a Myklestad method [6], 
a Ritz finite-element method [7], a modal approach based on a mixed varia- 
tional principle [8], a transmission matrix method [9], and a Galerkin 
finite-element method [10]. 
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In section 2 of this paper, as in [1] , tho governing equations of notion 
are written m two onto of aquations. Ona la a aat of nonlinaar ordinary 
diffarantlal aquations that govarns tha static equilibrium condition. Tha 
othas if a; sat of linear ordinary diffarantlal aquations with an unknown 
aiganvalua that govarns the dynamic behavior of snail perturbation notions 
Shout equilibrium* These diffarantlal equations are written In a mixed, 
spatial-derivative form, unlike tha dlaplacanent fomulation of [1]. 

Several differences between nixed and dlaplacanent formulations are then 
discussed in the context of the present analytical task. 

Then, In section 3, solution of the equations Is discussed. A differ- 
ent technique is used to solve for the nonlinear static equilibrium from 
that used to solve for the linear stability eigenvalues. The nonlinear 
static equllibriun Is solved by a collocation method for a mixed-order sys- 
tem of boundary value equations [11]. The software for this method Is 
available In a program called COLSYS [12]. For the linearised perturbation 
equations, several methods of solving for the eigenvalues are discussed. 

These include COLSYS [12], a generalised eigenvalue approach [13], a Block- 
Stodola technique [14], and one-step numerical integration techniques of 
various orders. For the present formulation the one-step methods, which 
appear to be the most promising, are used. 

Finally, In section 4, numerical results are presented for comparison 
with published data and for study of convergence properties associated with 
the one-step methods. Numerical results are also presented for stability 
eigenvalues of a hypothetical nonuniform blade. 

2. The Governing Equations 

The governing equations used herein were essentially derived In [15], 
in which the Houbolt and Brooks [16] equations are extended into the range 
of geometric nonlinear behavior. The integro-partial differential equa- 
tions are transformed into ordinary differential equations by first express- 
ing the vector of unknowns, z, as 

z(x,t) - z(x) + z(x,t) (1) 

Here 2 is the static equilibrium part and 2 is a small, perturbed part. 
The equations are linearized in z and converted to an eigenvalue problem 
with the transformation 

z'.x.t) - z(x)e Xt (2) 

where 

X - a + i<o (3) 

The static equilibrium equations (with the sign conventions and nomenclature 
defined in [15] are given in appendix A. From Equations (Al) through (A10), 
and (A17), the state vector, 2, for the nonlinear, static equilibrium 
behavior is written as 

l - [v y M z v l V z M y w l M x 4 ?J T (4) 

and Equations (A13) , (A14) , and (A15) define the aerodynamic loading for 
the hovering flight condition. 
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The linearised squat ions govern Ing small perturbation notions from the 
static equilibrium state defined by the solution for the equations in appen- 
dix A are given in appendix B. The state variables, z, for the perturbed 
state are defined in Equations (Bl) through (BIO), and (B14) through (B19). 
The vector a is 

a - |v y v C V # M y w a H x ♦ i> 
and the equations may be represented in matrix form as 


a a n 

V x u J 


(5) 


s' - (A + XB)i (6) 

In the derivation cf the equations in appendices A and B, terms of 0(e 2 ) 
have been dropped with respect to unity where e is the order of bending 
and torsion rotations, g, 0, and The ordering scheme outlined in [15] 
has been followed as closely as possible. As in [15], an exception has 
been made to include 0(e 2 ) terms in the torsion equations that are uncou- 
pled from the bending equations but are known to influence the uncoupled 
torsion frequency. It is not possible, however, to be completely consistent 
in ordering schemes in either a mixed formulation for the differential 
equations or in a displacement formulation. Furthermore, as will be shown 
in the next section. It is not always possible to reach complete agreement 
between the mixed and displacement formulations in the decision of which 
terns to retain, even when trying to ' follow the same guidelines on neglect- 
ing higher-order terns. 


An exception to the guidelines in [15] on the ordering of terms has 
been to retain complete trigonometric expressions involving the variable 4). 
In [15], cos(6 + +) and sln(6 + 4) were approximated by the appropriate 
expansions in terms of cos 6, sin 6, and neglecting terms of second and 
higher powers of 4 . In this paper it has been convenient to retain the 
complete trigonometric expression in all equations except the equations 
relating rotation derivatives and moments which are presented below. - No 
significant differences in the results of [1] and this paper will occur as 
a consequence of using this convenient trigonometric form. 


When the equations are formulated as a system of first-order differen- 
tial equations (i.e., a mixed formulation) several attractive features 
become apparent. These include the simplicity of the form of the governing 
differential equations, the absence of the derivatives of the elastic char- 
acteristics, the simplicity of applying numerical integration techniques, 
and of handling boundary and interface conditions. The mixed formulation 
was particularly convenient in this application since the integro- 
differential terms could easily be rewritten as differential equations 
[note Eq. (A17), (B14) , and (B19)]. 


Concerning retention of specific terms, the main difference between a 
mixed-order formulation and displacement formulations involves the equations 
relating derivatives of rotations to moments which are 


JC’J _ fajC .7 it j 

U’J L aic - ai6 


-aj + a^J 
-a* - 2ai$ 


az 

ax 


&*- 


«x 

MzJ 


+ 


(•aV z «> 


I cos(6 
sin(0 


+ $) 
+ 4 >) 


(7) 
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These aquations lead to Equations (A4) and (A8) in appendix A and Equa- 
tions (B4) and (B8) In appendix B. Reference [1] , which is a displacement 
formulation neglects the M contribution (i.e. , the single-underlined terms 
in Equation 7) In a consistint manner on the basis that torsional moments are 
at least one order of magnitude less than the bending moments. References [1] 
and [15] use a quasi-coordinate as the torsion variable as discussed in [17]. 
For that formulation, the torsion moments when used in the bending equations 
are written as Integrals of applied and inertial loads and are thus of 
0 (r 2 ). In the bending equations they are multiplied by quantities of 0(e) 
and are thus negligible. References [18] and [19] , which are also displace- 
ment formulations, use Lagrangian torsion variables and in a consistent 
manner retain the single-underlined terms. There, the torsion moments are 
written in terms of the first derivatives of $ and are 0(e) quantities. 

When they are multiplied by terms of 0(e) in the bending equations, they are 
not formally negligible and thus are retained. The single-underlined terms 
are also retained in the mixed fc relation presented here, although they could 
be neglected on the basis of the arguments in [1] . Results are presented in 
a later section which document the effect of the single-underlined terms. 

Another more important difference between mixed and displacement for- 
mulations is illustrated by the double-underlined terms of Equation (7). 

In [1]* [15], [18], and [19] the double-underlined terms were retained in 
deriving the fourth-order governing equations for the lead-lag and flap- 
bending behavior, but they were neglected in deriving the second-order 
torsion equation. These terms are consistently retained or neglected in 
the above references on the basis of essentially the same ordering scheme, 
and the retained terms result in a symmetric elastic stiffness matrix. For 
the present mixed formulation, however, there is no reason, a priori , to 
neglect these terms. We may expect these terms to have some observable 
effect, especially when there is significant static equilibrium deformation, 
such as when 8 is large. This effect will be Illustrated in the next 
section. 

Previously, we discussed some advantages of the mixed formulation. There 
are also some drawbacks to use of mixed formulations. The matrix operator for 
the structural terms in the equations no longer appears in symmetric form, and 
the number of variables in the state vector is several times that of displace- 
ment formulations. For some applications, however, the convergence is more 
rapid for mixed formulations than for displacement formulations [20 and 21], 
thereby offsetting somewhat the latter disadvantage. Further, force resultants 
are obtained with the same degree of accuracy as the displacements in the 
mixed formulation. The disadvantage of the nonsymmetrlc form of the structural 
operater can be somewhat offset by selection of solution techniques that do 
not make use of matrix symmetry or positive definiteness. Basically, the 
choices considered here were whether to use collocation methods, generalized 
eigensystem approaches, or one-step integration methods. 

3. Solution of the Systems of Equations 

For obtaining the nonlinear, static equilibrium solution, the colloca- 
tion method [11] was used since the software [12] was readily available 
and since the method ensures a high accuracy. The program is capable of 
treating general systems of nonlinear multipoint boundary value problems up 
to order four with a variety of options available to the user. COLSYS pro- 
duces an approximate solution to a user-specified accuracy using a polynomial 
spline which can be evaluated at any point within the domain of the problem. 

The Inertial, geometric, and elastic properties of the blade can be expressed 
as functions of the axial coordinate in this approach, and the high-order. 
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spline-fit method leads to almost as high a degree of accuracy at noncollo- 
catlon point a as at the collocation points. Further, trigonometric expres- 
sions involving state variables, such as sin (6 + $) , can be expressed 
exactly without the small angle assuaptions on $ that are often made. 

The collocation software package, COL SYS, was written to solve a set 
of nonlinear multipoint boundary value problems. As such, it can also be 
applied tc eigenvalue problems using the approach outlined in [22]. After 
investigation of this approach for the eigenvalue problem, it was discarded 
for various reasons. In particular, since the complex roots doubled the 
number of required governing equations, the existing code could no longer 
accommodate these cases without modification. A further limitation is that 
a priori estimates of the desired complex eigenvectors and eigenvalues 
must be provided. 

Next, the generalised eigenvalue approach was investigated. In this 
approach. Equation (6) is discretized using a finite-difference approxima- 
tion. The one chosen here involved the central difference approximation 
given by: 

- ‘l - (A + XW 1+l/2 + « 1+1 > (») 

where the subscripts denote the locations at which the variables are 
evaluated. 

This results In a generalized eigenvalue problem of the form 

• (A + XB)Z - 0 . (9) 

The dimensions of the matrices A and B are 16k x 16k, where 16 is the 
dimension of the matrices A and B in Equation (6), and k is the num- 
ber of segments in the discretization. In preliminary investigations the 
subroutine EIGZF [13] was used to obtain all the eigenvalues for this 
generalized eigenvalue problem. It was found that as many as 100 segments 
would be required to obtain 3-place accuracy for the second eigenvalue. 
Unfortunately, the lack of symmetry prevents the banded structure of the 
matrices from being exploited; thus, a storage problem is created. 

One alternate approach is to find selected eigenvalues. The Blcck- 
Stodola method [14] uses a block inverse-iteration algorithm to find the 
first few eigenvalues. In those cases which were tested for this method, 
no more than the first eigenvalue could be computed accurately. Here, 
because of the lack of symmetry, there is no underlying variational princi- 
ple that can be employed as in [14]. Although it is possible to symmetrize 
one of the two matrices and obtain the hypotheses in [14], this results in 
a significant degree of fill-in. As a consequence, this approach was not 
pursued. 

The one-step methods using transfer matrices similar to the Myklestad 
method seem to be the most promising. The term one-step method refers to 
those methods which depend only upon station i to obtain solutions at 
station 1 + 1. Higher-order terms of a Taylor expansion can be used with 
the one-step method to speed convergence. The method can be described as 
follows. A Taylor expansion of z ( x i+l ) about is 

z(x 1+1 ) - zfr^ + z'Cx^h + z"(x i )h 2 /2 + 2? ,, (x i )h 8 /6 + . . . (10) 
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Now, in Equation (6), lot 


h " x l+i " x i 


C - A + XB 


(ID 


( 12 ) 


A first-order method eon be obtained from the first two terns of Equa- 
tion (10), yielding 


1+1 


*i + 


M i*i 


(13) 


where Equation (6) la used to eliminate z'(xi). Equation (13) has 0(h 2 ) 
local truncation error. A second-order method can be similarly obtained 
from the first three terms in Equation (10) yielding 


* 1+l ■ « i + be i * 1 + (h 2 /2) (CJ* 1 + Cjjtj) 


(14) 


which has 0(h s ) local truncation error. Here Equation (6) is used to 
eliminate z'(x^) and z”(x±). Both Equations (13) and (14) may be rewritten 
as Z£tx - Tjz^, where is a suitably chosen matrix evaluated at . xi, 
so that 


z i+i " (I + hC i )z i 

z i+i " [1 + “i + ( h2 / 2 >( c i + 


(15) 

(16) 


The form of depends on the order of the method. Hence, z n ■ Tzo where 
T - T n -xT n _ 2 **»To. The homogenous boundary conditions at Xq and Xq are then 
used to reduce matrix T to a smaller (6 x 6) matrix f whose determinant 
det (?) is a polynomial In X. The desired eigenvalues correspond to those 
X's that produce zero determinant for T. In cur Implementation of the 
one-step methods, subroutine ZANLYT [13] Is used to find one or more of the 
complex roots of the real golynomial det(f). Subroutine LEQTIC [13] Is 
used to decompose matrix T into L and U factors. The product of the 
pivots produced from this factorization gives the determinant of T. Finally, 
whenever Is used, it Is approximated by the backward difference 

C 1 “ (C i ' C i-i )/h (17) 

and Cq is assumed to be zero. In any one-step method, care must be taken 
to* define the C' terms properly at those points where cusps or discon- 
tinuities in the axially varying properties occur. 

If there are axial variations in inertial, geometric, and elastic 
properties or In the tension, the term Cjz^ in Equation (16) can have sig- 
nificant impact on the speed of convergence. Rotor blades have all of 
these axial variations. Therefore, the method illustrated by Equation (16) 
was used In this paper. As shown in the next section the nuinber of seg- 
ments required to produce accurate results increases substantially if the 
Cl z i term is omitted from Equation (16). 
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4. Numerical Results and Discussion 


Numerical reeult* vara obtained using an IBM 360/67 computer and are 
presented in this section for several rotor-blade configurations. These 
results are intended to serve primarily three purposes: (1) to compare 

with published results and thus, at least partially, validate the present 
computer program; (2) to study the convergence properties of one-step 
methods; and (3) to present some new results for a nonuniform blade which 
may serve as a reference problem for future analytical studies. The 
numerical values of the various inertial, elastic, and geometric properties 
are presented in Table 1 for configurations to be compared with results 
from [1] and [9] and for the hypothetical nonuniform blade. 

4.1. Comparison with published results 

In this section numerical results are presented to compare with some 
published data. No attempt will be made to compare with all of the many 
available numerical results in the literature. Instead, we will focus on 
results from an in vacuo configuration in [9] and the aeroelastic stability 
results of [1]. Because of the very good agreement between the present 
results and those published, it was felt that a preeentatlon of the results 
in tabular form would facilitate comparison with published data and under- 
standing of the effects of the underlined terms in Equation (7) on the aero- 
elastic stability results of [1]. For serving the latter purpose we 
designate two constants kl and k2. For kl (or k2) « 1, the single- (or 
double-) underlined terms in Equation (7) are Included. For kl (or k2) * 0, 
the single-(or double-) underlined terms in Equation (7) are deleted. 

As pointed out above, consistent application of the ordering scheme in 
a displacement formulation, such as in [1], may lead to discarding terms 
that may not necessarily be negligible in an analogous ordering scheme for 
a mixed formulation as in thi" paper. In order to match results of [1], it 
was necessary to experiment with deleting and retaining the underlined terms 
in Equation (7). One example of the dilemmas faced in trying to be con- 
sistent is apparent when the torsion moment M x is written as an Integral 
of applied and inertial loads. In this case it is clearly 0(e 2 ), as noted 
in [15] and explained in [17]. However, when Mx is written in terms of 
the derivatives of rotations, it contains one term, which is 0(c). If Mx 
is regarded as 0(e), the single-underlined terms must remain (kl - 1); if 
M x is regarded as 0(e 2 ), these terms should be neglected for the sake of 
consistency (kl * 0) . A second example Involves the double-underlined 
terms. When the moment components My and M 2 appear in the bending equa- 
tions in a displacement formulation, they appear multiplied by 0(1) 
quantities. On the other hand, in the torsion equation they are multiplied 
by 0(e) quantities. Because terms up through 0(e 2 ) are retained in all 
equations in [1], it is consistent to take only the dominant terms in 
expressions for My and M z in the torsion equation but retain the more 
accurate expression, including double-underlined terms, for My and ^ in 
the bending equations. In the present mixed formulation, however, dropping 
of the double-underlined terms (k2 ■ 0) in Equation (7) would result in 
matching the torsion equation of a displacement formulation equation but 
result in an oversimplified set of bending equations. Retention of the 
double-underlined term (k2 ■ 1) results in matching the bending equations 
of a displacement formulation, but adds additional 0(e s ) terms that are 
inconsistent in the torsion equation of the displacement formulation. Thus, 
a consistent set of equations in a mixed formulation may, when transformed 
to a displacement form, result in neglect or retention of terms that would 
be differently treated in a set of equations written consistently in 
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TABLE 1. VALUES Of INERTIAL, ELASTIC, AMU GEOMETRIC PARAMETERS FOR 
COMPARISON OF RESULTS WITH (1) AMD [9]. 


(a) Reference [1] Configuration 
p„cR/a » 5/(6*) 
bc/(*R) - 0.1 


c/R ■ ir/40 

a » 2w 

c^/a - 0.01 

Ely/CaOV) - 0.014605 <u> w - 1.150) 

( 0.026787 (u» * 0.70, soft in-plane case) 
0.16696 (bv * 1.50, stiff in-plane case) 
GJ/(nfl*R % ) - 0.0056732 («^ - 50) 

W R " *a 2 /* - 0*025 
fc^i - a A - a - 0 
(k A /k,) 1 - 1.5 
6 » 0, 0.3 rad 

(b) Reference [9] Configuration 

0*0 0 - ir/4 rad 

P» - 0 GJ « 9000 lb in. 2 


R * 40 in. k 2 ^ ■ k 2 ^ ■ 1 in. 2 

m - 0.000125 lb sac 2 in. -2 a. - k. - 0 

A A 

EI y - 25000 lb in. 2 e - /l in. 

Elg - 7500 lb in. 2 

(c) Tapered Case 

All properties at blade root identical to those listed above in part (a) 
for the soft-in-plane case. All properties are constant except Ely, 
EI Z , GJ, and a, which are multiplied by the factor (1 - O.lx). 
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d lap lac— at form. We expect differences such as that* to be nail. 
Otherwise, we mac eoneluda that whatever has baaa taken aa an ordaring 
achaaa auat itself ba somehow Inconaiatant. 

A aaapla of raaulta ganaratad for tha simplified configuration from 
[9] ia praaantad in Tabla 2. Tha rotor apaad la taro in thla caaa, and 
tharafora thara la no aerodynamic loading and no static equilibria defor- 
aatlon. Thara ia navarthalasa, an offset between tha blade aass canter and 
shear canter and nonaaro pitch angle (tha properties are given in Tabla 1). 
Thus, all degrees of freedom are coupled. Tha number of segments used by 
tha example from [9] is unknown. However, it is clear that tha present 
results baaed on a second-order, one-step method (O' • 0 for this case, 
since tha bean ia uniform with no tension) are tending towards those of 
[9] as tha mnriber of segments is increased. With only 24 segments, the 
first mode is within 0.12 of Murthy's result. 

TABLE 2. COMPARISON WITH RESULTS FROM 
[9] FOR FREE VIBRATION FREQUFNCIES OF 
BLADE WITH MASS CENTER OFFSET. 


u> (rad/ sec) 


Segments 

Mode 1 

Mode 2 

Mode 3 

16 

30.7552 

53.6968 

179.8088 

24 

30.7962 

53.7691 

182.3980 

32 

30.8107 

53.7947 

183.3519 

40 

30.8174 

53.8065 

183.8023 

Ref. [9] 

30.8295 

53.8277 

184.6175 


When the rotor speed is nonsero, the presence of steady aerodynamic 
loads and, in some cases. Inertial loads, causes significant static defor- 
mation which must be taken into account properly in order to obtain correct 
stability eigenvalues [1]. Results from the present analysis (COL SYS) are 
presented in Table 3 along with those of [1] for both roft-in-plane (u>v ■ 0.7) 
and stiff-in-plane (ul. ■ l.S) cases with_ 0-0 and 0.3. The quantities 
tabulated are tip deflections of v, v, ?, and the accuracy is specified to 
be four significant figures. The solution is neither difficult to obtain 
nor particularly time consuming. The agreement is quite good regardless of 
the choice of kl and k2. We conclude from this that the static equilibrium 
solution is not strongly affected by the presence of the underlined terms 
in Equation (7) for the limited range of values presented here. Such is 
not the case fox the stability eigenvalues, however. 

Table 4 presents the stability eigenvalues based on a second-order, 
one-step method with C* terms included and 16 segments (converged to 
3 significant figures). Again, the agreement is very good, especially at 
8-0. The presence of the single-underlined terms (kl - 1) has no effect 
on the results for 0-0 and very little effect for 6 - 0.3. The double- 
underlined terms (k2 ■ 1) raise the torsion frequency at 6 - 0.3 from its 
value at 0-0 Instead of lowering it as indicated in [1]. The effects 
of k2 on lead-lag damping ere minor. The difference in the torsion fre- 
quency due to the k2 terms is about 22. In a displacement formulation, 
these terms of second degree in bending curvatures are neglected in the 
torsion stiffness (i.e., third-degree terms, such as *v"*, 4v"w", etc., 
are neglected in the second-degree torsion equation). This is consistent 
in displacement formulations such as {1], [15], [18], and [19]. It is not 
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TABU 3. VALOIS OF STATIC IQUILmiOM 
DISPLACEMENTS AT BLADE TIP. (Id - 0, 

* k2 ■ 1 for COLTS reauAta) 


v/R 

9/R 


(a) • 0.70, 8 ■ 

0.0 

COLSTS -0.002321 

0.0 

0.0 

Raf. [1] -0.002326 

0.0 

0.0 

(b) Ny ■ 0.70, 8 ■ 

0.3 

COLSTS -0.03943 

0.09315 

-0.01332 

Raf. [1] -0.03940 

0.09314 

-0.01336 

(c) « v - 1 

.50, 8 - 

0.0 

COLYST -0.000522 

0.0 

0.0 

Raf. [1] -0.000522 

0.0 

0.0 

(d) - l 

.50, 8 - 

0.3 

COLSTS -0.03139 

0.09432 

-0.01164 

Raf. [1] -0.03129 

0.09412 

-0.01214 


TABU 4. COMPARISON OF STABILITY EIGENVALUE RESULTS OBTAINED BY 
USING VARIOUS VALUES OF kl and k2 (16 aegaenta). Eigenvalue a are 
given per unit Q. 


(a) 8 • 0 

kl 

k2 

Lead-lag 

Flapping 

Tore Ion 

My - 

0.7Q 







0 

0 

-0.0010 

tO. 69291 

-0.3235 

±1.07881 

-0.3615 

±4.98091 

0 

1 

-0.0010 

±0.69291 

-0.3235 

±1.07881 

-0.3615 

±4.98091 

1 

1 

-0.0010 

±0.69291 

-0.3235 

±1.07881 

-0.3615 

±4.98091 

Raf. 

Cl] 

-0.0011 

±0.70141 

-0.3245 

±1.07511 

-0.3622 

±4.98751 

“\T - 

l.SQ 







0 

0 

-0.0011 

±1.49381 

-0.3235 

±1. 0767i 

-0.3617 

±4.98221 

0 

1 

-0.0011 

±1.49361 

-0.3235 

±1.07671 

-0.3617 

±4.98221 

1 

1 

-0.0011 

±1.49381 

-0.3235 

±1.07671 

-0.3617 

±4.98221 

Raf. 

[1] 

-0.0011 

±1.50021 

-0.3246 

±1.07411 

-0.3625 

±4.98881 





(b) 8 - 0. 

.3 



<*v - 

0.70 







0 

0 

-0.0230 

±0.68231 

-0.3145 

±1 .07491 

-0.3552 

±4.95921 

0 

1 

-0.0228 

±0.68191 

-0.3146 

±1.07461 

-0.3564 

±4.97201 

1 

1 

-0.0226 

±0.68151 

-0.3147 

±1.07431 

-0.3574 

±4.98361 

laf. 

(11 

-0.0249 

±0.69311 

-0.3117 

±1.06761 

-0.3628 

±4.96371 

mv - 

1.50 







0 

0 

-0.0641 

±1.61431 

-0.2821 

±0.95471 

-0.3520 

±4.95501 

0 

1 

-0.0665 

±1.56981 

-0.2792 

±0.95451 

-0.3494 

±5.04711 

1 

1 

-0.0667 

±1.56631 

-0.2789 

±0.95451 

-0.3491 

±5.05241 

Raf. 

m 

-0.0632 

±1.58231 

-0.2821 

±0.95071 

-0.3510 

±4.91501 
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obvious, hove v fir , chat these third-degree terns are the only ones that have 
an observable effect on the trends. Ths only way to resolve that question 
is to examine results froa a consistent analysis in which O(e’) terms are 
neglected with respect to unity, which has not yet been developed. 

4.2. Convergence of one-step methods 

He now address the subject of convergence for the one-step methods 
including the C* terms. In Table 5, results obtained for the soft-in-plane 
case at 0 ■ 0.3 are tabulated as a function of the number of segments. 

The convergence is fairly rapid at first and tapers off as the converged 
value is approached. The smallest eigenvalues are generally within 1% 
for 16 segments. 


TABLE 5. CONVERGENCE OF THE ONE-STEP METHOD BASED ON EQUATION (16) 
FOR VARIOUS NUMBERS OF SEGMENTS <% - 0.70, 0 - 0.3, kl - 0, k2 - 1) 
Eigenvalues are given per unit 0. 


Segments 

Lead-lag 

Flapping 

Torsion 

4 

-0.0159 ±0.7394i 

-0.3099 

±1.20251 

-0.3365 ±4. 89601 

8 

-0.0234 ±0.6902i 

-0.3108 

±1.09291 

-0.3523 ±4. 95481 

12 

-0.0230 ±0.68401 

-0.3136 

±1. 07911 

-0.3553 ±4.96761 

16 

-0.0228 ±0.68191 

-0.3146 

±1.0746i 

-0.3564 ±4.9720i 

20 

-0.0227 ±0.68101 

-0.3151 

±1.07251 

-0.3569 ±4.9739i 

24 

-0.0227 ±0.68051 

-0.3153 

±1. 07151 

-0.3572 ±4. 97501 


In the one-step method, the C' terms may be neglected in some appli- 
cations. This is certain to have an adverse effect on convergence for a 
rotating beam, however. Even for a beam with uniform properties, the per- 
turbation equations will have variable coefficients due to tension and 
static equilibrium terms. Thus, it is important to study the convergence 
of the stability eigenvalues as a function of the number of segments for 
methods which do not use the C' terms. In Table 6, the lead-lag and flap 
eigenvalues are presented for three methods versus the number of segments. 
Method (1) is the complete second-order method with the C* terms included 
(Eq. (16)). Method (2) is the second order method without the C' terms. 
Method (3) is the first-order method (Eq. (15)). There is little differ- 
ence in the rate of convergence for Methods (2) _ad (3). The inclusion of 
the C' terms, however, results in about an order of magnitude reduction 
in the number of segments required for convergence. Hundreds of segments 
may be needed for convergence when C 1 is neglected [23]. It should be 
noted that [6] uses only the tension terms in evaluation of C' . While 
that should result in a substantial improvement in convergence over Meth- 
ods (2) and (3), it would still be inferior to Method (1) if the beam has 
portions of even modest nonuniformities in stiffness or inertial properties. 

4.3. Results for a nonuniform blade 

We have included results for a hypothetical nonuniform blade in Table 7. 
All properties are the same as those used in the configuration for compari- 
son with [1], except that Ely, EI Z , GJ, and m are equal to their values 
for the uniform blade times the factor (1 - O.lx). Convergence for the 
nonuniform blade is a little slower, and all frequencies are increased from 
the values obtained the uniform blade. The lead- lag damping is only 
slightly decreased by the presence of taper for this case. This small 
effect is not surprising, since fundamental bending- torsion coupling 
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TABLE 6. EFFECT OF C]zi TERMS ON CONVERGENCE OF LABILITY EIGEN 
VALDES USING ONE-STEP METHODS. <6 - 0.3, a* - 0.70, kl - 0, 
k2 ■ 1). Eigenvalues are given per unit 0. 


Segments Method (1) Method (2) Method (3) 


Lead- lag 


8 

-0.0234 ±0.69021 

-0.0201 

±0.76991 

-0.0181 ±0.83981 

16 

-0.0228 ±0.68191 

-0.0216 

±0.73081 

-0.0203 ±0.7500i 

24 

-0.0227 ±0.68051 

-0.0220 

±0.71521 

-0.0211 ±0.72441 

32 


-0.0221 

±0.70681 

-0.0214 ±0.71241 

48 


-0.0223 

±0.69801 

-0.0218 ±0.70091 

64 


-0.0224 

±0.69351 

-0.0220 ±0.69531 

80 


-0.0224 

±0.69081 

-0.0222 ±0.69201 

96 


-0.0225 

±0.68891 

-0.0222 ±0.68991 


Flap 


8 

-0.3108 

±1.09291 

-0.2849 

±1.13671 

-0.2922 

±1. 14311 

16 

-0.3146 

±1.07461 

-0.3008 

±1.10601 

-0.3034 

±1.10211 

24 

-0.3153 

±1.07151 

-0.3059 

±1.09441 

-0.3075 

±1. 09031 

32 



-0.3085 

±1.08841 

-0.3095 

±1.0847 i 

48 



-0.3110 

±1.08211 

-0.3116 

±1.07931 

64 



-0.3122 

±1.07891 

-0.3127 

±1.07671 

80 



-0.3129 

±1.0770i 

-0.3133 

±1.0751 i 

96 



-0.3134 

±1.0757 i 

-0.3137 

±1.07471 


TABLE 7. 

STABILITY EIGENVALUES FOR A NONUNIFORM BLADE (kl » k2 - 1 
6 “ 0.3). Eigenvalues are given per unit 0. 

Segments 

Lead-lag 

Flap 

Torsion 

16 

-0.0230 ±0.70011 

-0.3138 ±1. 08191 

-0.3564 ±5.07981 

24 

-0.0229 ±0.69871 

-0.3144 ±1. 07901 

-0.3567 ±5.08371 

32 

-0.0228 ±0.69831 

-0.3146 ±1.07811 

-0.3564 ±5.0851 i 


parameters [1] are not changed significantly. No significant computational 
penalty in the calculation of the nonlinear static equilibrium and the 
stability eigenvalues is incurred because of nonuniformity. 

5. Concluding Remarks 

Nonlinear equations for static equilibrium deformation and linearized 
perturbation equations for small motions about equilibrium, from which 
stability eigenvalues can be obtained, are written in a first-order, spatial 
derivative form. These equations are solved by COLSYS [12] and one-step 
Integration techniques, respectively. Numerical results for uniform blades 
obtained from the equations are compared with published results [1] and [9]. 
and used to study the convergence of the one-step methods. Results are 
also presented for a hypothetical nonuniform blade. In the course of this 
study several conclusions have emerged. 
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1. The equations of the present mixed formulation are general enough 
to treat nonuniform pretwisted rotor blades with chordvise offsets between 
elastic center, mss center, and shear center. Certain higher order cross- 
section integrals are neglected, and symmetry about the major cross-section 
axis is assisted [15]. Neither the equations nor the solution methods, in 
their present form, apply to the forward-flight problem. 

2. COL SYS is used to solve the nonlinear static equilibrium equations 
with the present mixed formulation. The calculation of the static equilib- 
ritsa solution is neither particularly difficult nor time consuming even 
though nonlinearities and nonuniformities are involved. Thus, there is no 
need to limit the static equilibrium solution to some linear or otherwise 
ad-hoc estimate. 

3. There are slight differences in the numerical results obtained from 
mixed and displacement formulations of rotating blades with geometric non- 
linearity. The differences stem from the equations that relate moments to 
derivatives of bending rotations. A consistent ordering scheme applied to 

a mixed formulation may not produce exactly equivalent equations when 
applied to a displacement formulation. The main difference in the results 
is in the magnitude of the torsion frequency. The difference is of the 
order of rotations squared with respect to unity (about 2Z). The present 
analysis yields the result that torsion frequency increases with increased 
pitch angle in direct contrast to results in [1]. 

4. One-step, second-order integration methods appear to be a viable 
means of calculating the stability eigenvalues. In order to obtain good 
convergence it is necessary to include the entire second-order term, for 
rotor blades, for which tension force, elastic, and inertial properties all 
may vary along the length of the blade. 

5. The terms involving the torsion moment ^ in Equation (7) were 
found to be negligibly small for the limited range of parameters investi- 
gated. The double-underlined terms do affect a basic trend, as described 
in conclusion (3) . This may indicate the need for consistently incorporat- 
ing terms of the next higher order in elastic rotations. Unfortunately, 
the only way to ascertain the correctness of this assertion is to compare 
results from such an analysis, which does not yet exist, to those of this 
paper or [1]. 


13-13 


-I ■ *> *>. *« » *}i**-*«. am.**- ur .*.**.»•* 


REFERENCES 



s 


l 

* 

l 


i 

5 

1 

i 

i 

\ 


I 

? 

t 

t 

f 


\ 

t 


1. Dewey H. Bodges and Robert A. Ormiston, Stability of elastic bending 
and torsion of uniform cantilever rotor blades in hover with variable 
structural coupling, NASA IN D-8192 (1976). 

2. Robert A. Ormiston and Dewey H. Hodges, Linear flap-lag dynamics of 
hingeless helicopter rotor blades in hover, J. American Helicopter Soc . 
i2. No. 2, pp. 2-14 (April 1972). 

3. Dewey H. Hodges and Robert A. Ormiston, Nonlinear equations for bend- 
ing of rotating beams with application to linear flap-lag stability 
of hingeless rotors, NASA TM X-2770 (1973). 

4. William F. White, Jr. , Importance of helicopter dynamics to the mathe- 
matical model of the helicopter, AGARD Flight Mechanics Panel Special- 
ist Meeting, Langley Research Center, Paper 22 (November 1974). 

5. William F. White, Jr., and Raymond E. Malatlno, A numerical method for 
detexminlng the natural vibration characteristics of rotating nonuni- 
form cantilever blades, NASA TM X-72,751 (1975). 

6. J. R. Van Gaasbeek, T. T. McLarty, and S. G. Sadler, Rotorcraft flight 
simulation, computer program C-81, Volume 1 — Engineer's Manual, 
USARTL-TR-77-54A (October 1979). 

7. Dewey H. Hodges, Vibration and response of nonuniform rotating beams 
with discontinuities, J. American Helicopter Soc . 24, No. 5, 

pp. 43-50 (October 1979) . 

8. K-W. Lang and S. Nemat-Nasser, An approach for estimating vibration 
characteristics of nonuniform rotor blades, AIAA J . 17, No. 9, 

pp. 995-1002 (September 1979) . 

9. V. R. Murthy, Dynamic characteristics of rotor blades, J. Sound .and 
Vibration 49, No. 4, pp. 483-500 (1976) . 

10. P. P. Friedmann and F. Straub, Application of the finite-element 
method to rotary-wing aeroelasticity, J. American Helicopter Soc . 25, 
No. 1 (January 1980). 

11. U. Ascher, Jay Christensen, and R. 0. Russell, A collocation solver 
for mixed order systems of boundary value problems, TR 77-13, Dept, of 
Computer Science, University of British Columbia, Vancouver, Canada 
(1977) (to appear in J. Math. Comp .). 

12. U. Ascher, Jay Christensen, and R. 0. Russell, C0LSYS-A collocation 
code for boundary value problems, Proc. of Conf. for Codes for - 
Boundary Value Problems in Ordinary Differential Equations, Houston, 
Texas (1978) . 

13. The International Mathematical and Statistical Library (IMSL) , Houston, 
Texas (1975). 

14. Stanley B. Dong, A Block-Stodola eigensolution technique for large 
algebraic systems with non-symmetrical matrices. International J. for 
Numerical Methods in Engineering 11, pp. 247-267 (1977). 


13-14 



16. John C. Houbolt and George W. Brooks* Differential equations of motion 
for combined flapvlse bending* chordwlse bending and torsion of twisted 
nonuniform rotor blades, NACA Rep. 1346 (1958). 


17. Dewey H. Hodges* Robert A. Orals ton* and David A. Peters, On the non- 
linear deformation geometry of Euler -Bernoulli beams, NASA TP-1566 
(1980). 

18. Krishna Rao V. Kaza and Raymond G. Kvatemik, Nonlinear aeroelastic 
equations for combined flapwise bending, chordwlse bending, torsion, 
and extension of twisted nonuniform rotor blades in forward flight, 
NASA TM-74059 (1977). 



19. A. Rosen and P. Friedmann, Nonlinear equations of equilibrium for 
elastic helicopter or wind tunnel blades undergoing moderate deforma- 
tion, UCLA-ENG-7718 (Revised edition) (June 1977). 

20. Ahmed K. Noor and Wendell B. Stephens, Mixed finite-difference scheme 
for free vibration analysis of noncircular cylinders, NASA TN D-7107 
(1973). 

21. Ahmed K. Noor and Wendell B. Stephens, Comparison of finite-difference 
schemes for uialysis of shell of revolution, NASA TN D-7337 (1973). 

22. Herbert B. Keller, Numerical solution of two-point boundary value 

problems, region Conf . Series in Applied Mathematics, SIAM J. Numeri- 
cal Analysis 24 (1976). ‘ 

23. T. J. McDaniel and V. R. Murthy, Bounds on the dynamic characteristics 
of rotating beams, AIAA J . 15, No. 3, pp. 439-442 (March 1977). 

24. Dewey H. Hodges, Torsion of pretwisted beams due to axial loading, 

ASME J. Applied Mechanics 47, No. 2, pp. 393-397 (June 1980). 


13-15 


Steady State Equation* 


The static equilibrium equations are given in the following first- 
order fora with the sign conventions end nomenclature of [IS]. As in [15], 
section properties of higher order, B t * and B 2 *. warping rigidity and shear 
center offset terms, and Cx and Ci* are neglected. In presenting these 
equations, care has beeh taken to neglect terms that are 0(e 2 ) compared 
to unity. To avoid confusion with d/dx, the primes have been removed from 


the subscripts of Ely, EI S , and V x . 


Lead- lag equations: 


Vy - -afl 2 [v + e cos(6 + $)] - Ly 

(Al) 

- -V y + V x ; + mefl 2 * cos (8 + ♦) 

(A2) 

v* • 5 

(A3) 

c* - (a x c - a^M* + (-a x + a % $)My + (a 2 - 2a l ?)M r + 

CO8(0 + $)/Z 0 
(A4) 

Flap equations: 


% " "4 

(A5) 

- V z - V X B - men 2 '-: sin (6 + $) 

' (A6) 

w* - 6 

(A7) 


V - (a,C - + (-a, - 2a x *)My + (a x - aj)^ + e A V x sln(9 + ♦)/Z 0 

(A8) 


Torsion equations: 

M* • -V 2 C + VyB + mefl 2 v sin (8 + $) + mQ 2 (k£ - lc^JcosCe + $)sin(8 + 5) - M 

(A9) 

- 1/GofMjj + My(C + a e 0' sin P) + M Z (I - a s 0’ cos 6) - a 7 0’V x ] (A10) 
where, for brevity, 

G 0 - GJ + a 5 V x - a* M y sin 0 + a 6 cos 0 

Y 0 - Ely 

z 0 - EI Z - EA e A 2 (All) 

a x - -[ (1/Y 0 ) - (1/Z 0 )]sin 0 cos 0 
a 2 - (sin 2 0/Y o ) + (cos 2 0/Z Q ) 



a, - (coa* 6/Y 0 ) + (ain* 6/2 # ) 

*a 

*0 “ k A* B V z o 

•• " k A 2EA V Z 0 

a 7 - (EI z /Z 0 ) (k A * - EJ/RA) 

The term Involving a 7 In Equation (410) has been altered to agree with 
the recent reaulta of [24]. 


(All 

concluded) 


In the above equations, the tension In the In extensible blade Is given 
by the uncoupled expression 



The aerodynamic loading is derived In [1]. 
Is 


dx 

The steady component 


(A12) 


Ly « (p«»ac/2) [v^ 2 - O 2 x 2 (cd 0 /a) - Oxv^e + 9) ] (A13) 

L* - (p„ac/2) + fi 2 x 2 (0 + ♦+*)- n 2 xvg + (fl 2 xc8/2)] (A14) 

M . - 0 (A15) 

9 

The term 9 In Equation (A15) Is defined 


* 


/’ 


CB ' dx 


(A16) 


and Is represented In the analysis as the following additional ordinary 
differential equation 

*' - 53' (A17) 


with the initial boundary value specified by the values of £ and 0' at 
x(0). The Induced velocity is given by 


v A - agn(9 + 9 0 )flR(iro/8) [A + (12/iro) | 0 + + ,| - 1] (A18) 


where 

♦ 0 - ♦ at x/R - 0.75 (A19) 

and 


o ■ bc/irR 


(A20) 
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APPENDIX B 


Linearized Perturbation Equations 

The perturbation aquations given In [IS] have been converted to the 
linear ordinary differential aquations In first-order spatial derivative 
font with the eigenvalue X. Although not necessary for the method of 
solution finally adopted, they ere presented linear In X as well. 

Lead-lag equations: 

Vy • «{Xv* - eX$* sin 6 - ft 2 [v - e$ sin (6 + $)] 

+ 2ft [Xu - e(Xc cos 6 + X0 sin 6)]} - L y (Bl) 

n; - -Vy + V x i + CV X - «eft 2 x$ 8in(e + *) + 2meftXv cos e (B2) 

v’ - c (B3) 

V - (a x C - a*!)#* + (-a 2 + a^)My + (a 2 - + M x (a 1 C - a 2 3) 

+ [a % My - 2&Ji t - (e A /Z o )sln(0 + $>V X ]$ + e^ x cos(0 + $)/Z 0 (B4) 

Flap equations: 

V x » tn(Xw* + eX$* cos 0) - ^ (B5> 

% “ " V x 6 - BV X - meft 2 x$ cos(0 + ♦) - 2meftXv sin 0 (B6) 

w' - 8 (B7) 

B' - (a,t - a x 0)fi - (a s + 2a^)M w + (aj - a 4 *)M + M(a 3 C - a 3 0) 
x y z x 

- [2a x My + a 4 fi t - (e A /Z 0 )V x cos(0 +*)]$+ e k <i x sin(e + $)/Z 0 

(B8) 

Torsion equations: 

- ;V Z + V y B + BVy - me[Xv* sin(0 + *) - ft 2 v sin(0 + $) 

- fl 2 v cos (6 + $>$ - Xw* cos (ft + $) + 2ftXu sin 0] + mk m 2 X$* 

+ mft 2 (ki a - k^)* cos20 - (B9) 

(B9) 

- (1/GoH^ + [t + a e (0* + ?’)sin 0]dy + Myi; 

+ [6 - e,(0’ + $')cos 0]M Z + M Z B} (BIO) 

where the constants a^ and G 0 are defined in appendix A. 
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Equations (Bl) through (BIO) have been linearized in terms of X by 
using the relations 

v* ■ Xv (B 11) 

w* - Xw (B12) 

♦* - X$ (B13) 

Additional differential equations needed to supplement Equations (Bl) 
through (BIO) are 

- H' + Cfj' (B14) 

v*» - Xv* (B15) 

w*' - Xw* (B16) 

♦*' “ (B17) 

The perturbed tension and displacement in the axial direction are 


V’ - -2mflXv 
x 


(B18) 


u' - -H - eg - (e A /Z 0 )(My sin 9 


- M cos 6) (B19) 

z 

The aerodynamic load contributions [1] are 

Ly - (p s „ac/2){-nxv i $ - [2nx(c dfl /a) + (0 + ^v^Xv 

+ [2v i - Ox(0 + $)]Xw} (B20) 

L w ■ (p <x> ac/2){fi 2 x z ($ + <£) - fl 2 x(vg + 8v) + (Cl 2 xcB/2) 

+ [2tlx(6 + $) - v ]Xv - nxXw + (3/4)cftxX$ - (c/4)Xw*> 

(B21) 

- -(p w ac/2)[(c 2 /8)fixX$] (B22) 

Thus, Equations (Bl) through (BIO), (B14) through (B19) form the complete 
set of 16 coupled ordinary first-order differential equations required to 
calculate the stability eigenvalues of a nonuniform, pretvisted blade. 
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